Exercise 1: Invert a matrix


In [1]:
def matrix_is_square(M):
    if len(M) == len(M[0]):
        return True
    else:
        return False

class MatrixIsNotSquareException(Exception):
    pass

class MatrixIsSingularException(Exception):
    pass

def matrix_id(d):
    '''
    This function takes an integer d and
    returns the d*d identity matrix
    
    matrix_id(3)
    [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    
    '''
    return [i*[0]+[1]+[0]*(d-i-1)
            for i in range(d)]

def matrix_row_swap(M, upper_row, lower_row):
    '''
    >>> matrix_row_swap([[1,2],[3,4]],0,1)
    [[3, 4], [1, 2]]
    '''
    M[upper_row], M[lower_row] = M[lower_row], M[upper_row]
    return M
    
def matrix_row_multiplication(M, row_number, scalar):
    '''
    matrix_row_multiplication([[1,2],[3,4]],0,3)
    [[3, 6], [3, 4]]
    '''
    M[row_number] = [scalar * e for e in M[row_number]]
    return M

def matrix_row_addition(M, changed_row, mulitplied_row, scalar):
    '''
    matrix_row_addition([[1,2],[3,4]], 1, 0, 2)
    [[7, 10], [3, 4]]
    '''
    M[changed_row] = [sum(x) for x in zip(M[changed_row], 
                      [scalar * e for e in M[mulitplied_row]])]
    
    

def matrix_inverse(A):
    '''
    This function calculates the inverse B of a given matirx A
    such that AB = BA = I, where I is the identiy matirx of
    same dimension as A or B.
    
    >>> matrix_inverse([[2,1],[0,2]])
    [[1.0, -0.5], [0.0, 0.5]]
    
    >>> matrix_inverse([[1, 1, 0], \
                        [0, 1, 1], \
                        [2, 0, 2]])
    [[0.5, -0.5, 0.25], [0.5, 0.5, -0.25], [-0.5, 0.5, 0.25]]
    
    >>> matrix_inverse([[1, 1, 0, 1], \
                        [0, 1, 1, 1], \
                        [1, 0, 1, 1], \
                        [1, 1, 1, 1]])
    [[0.0, -1.0, 0.0, 1.0], [0.0, 0.0, -1.0, 1.0], [-1.0, 0.0, 0.0, 1.0], [1.0, 1.0, 1.0, -2.0]]
    '''
    if not matrix_is_square(A):
        raise MatrixIsNotSquareException
        
    B = matrix_id(len(A))
    
    # Make all elemens under the diagonal 0
    # For each column in the matrix
    for col in range(len(A)):
        # Swap the rows such that elements on the diagonal are not zero.
        if A[col][col] == 0:
            for row in range(col+1, len(A)):
                if A[row][col] != 0:
                    matrix_row_swap(A, col, row)
                    matrix_row_swap(B, col, row)
                    break
                else:
                    raise MatrixIsSingularException
        
        for row in range(col+1, len(A)):
            if A[row][col] != 0:
                scalar = -(A[row][col]/A[col][col])
                matrix_row_addition(A, row, col, scalar)
                matrix_row_addition(B, row, col, scalar)
    
    
    # Make all elemens over the diagonal 0 and the diagonal 1
    for col in range(len(A)-1, 0, -1):
        if A[col][col] != 1:
            scalar = 1/A[col][col]
            matrix_row_multiplication(A, col, scalar)
            matrix_row_multiplication(B, col, scalar)
            
        for row in range(col-1, -1, -1):
            if A[row][col] != 0:
                scalar = -(A[row][col]/A[col][col])
                matrix_row_addition(A, row, col, scalar)
                matrix_row_addition(B, row, col, scalar)
    
    
    return B   
    
if __name__ == '__main__':
    import doctest
    doctest.testmod()

Exercise 2: Visualize run-time


In [2]:
def fact2(n):
    '''
    This function takes a positive integer n
    and returns its factorial
    
    >>> fact2(0)
    1
    
    >>> fact2(1)
    1
    
    >>> fact2(8)
    40320
    
    '''
    
    assert type(n) == int, \
        'd should be an integer it is an '+type(n)
    assert n >= 0 , \
        'd should be a positive integer it is '+str(n)
    
    if n == 0: 
        return 1
    result = 1
    for x in range(1,n+1):
        result *= x
        
    return result


def prod(l):
    '''
    This function takes a list of integers and 
    returns the product of all elements
    
    >>> prod([1,2,3])
    6
    
    >>> prod([-1,2,-3])
    6
    
    >>> prod([3,0,2])
    0
    
    >>> prod([-3,5.5,1,7])
    -115.5
    
    '''

    assert type(l) in [list, range], \
        'l should be a list or range of numbers'
    
    if len(l) == 1:
        return l[0]
    if len(l) == 2:
        return l[0]*l[1]
    else:
        mid = len(l)//2
        return prod(l[:mid]) * prod(l[mid:])

    
def fact3(n):
    '''
    This function takes a positive integer n
    and returns its factorial
    
    >>> fact3(0)
    1
    
    >>> fact3(1)
    1
    
    >>> fact3(8)
    40320
    
    '''

    assert type(n) == int, \
        'd should be an integer it is an '+type(n)
    assert n >= 0 , \
        'd should be a positive integer it is '+str(n)

    if n == 0 or n == 1: 
        return 1    

    return prod(range(1,n+1))

        
if __name__ == '__main__':
    import doctest
    doctest.testmod()

In [3]:
import time
from functools import wraps

def timeit(func, calls=1):
    '''
    Decorator to time the execution time of a function
    '''
    @wraps(func)
    def timed(*args):
        measurements = []
        for call in range(args[-1]):
            start = time.perf_counter()
            result = func(*args[:-1])
            total = time.perf_counter() - start       
            measurements.append(total)
        avg_time = sum(measurements)/calls
        return result, total        
    return timed

In [4]:
import math
timed_fact2 = timeit(fact2)
timed_fact3 = timeit(fact3)
timed_math_fact = timeit(math.factorial)

x = []
time_fact2 = []
time_fact3 = []
time_math_fact = []

for exp in range(1,21):
    x.append(2**exp)
    time_fact2.append(timed_fact2(2**exp, 10)[1])
    time_fact3.append(timed_fact3(2**exp, 10)[1])
    time_math_fact.append(timed_math_fact(2**exp, 10)[1])

In [5]:
timed_math_fact(2**20, 10)[1]


Out[5]:
14.373267654998926

In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.loglog(x,time_fact2,':', label='fact2', )
plt.loglog(x,time_fact3, '--', label='fact3', )
plt.loglog(x,time_math_fact, '-.', label='math.factorial', )
plt.legend(loc=2)
plt.xlabel('Function input')
plt.ylabel('Execution time [s]')
fig = plt.gcf()
fig.set_size_inches(8,5)



In [14]:
skeleton = '{:^10.5} {:^10.5} {:^10.5}'
print(skeleton.format('frac2', 'frac3', 'mfrac')) 
for x in zip(time_fact2, time_fact3, time_math_fact):
    print(skeleton.format(*x))


  frac2      frac3      mfrac   
2.231e-06  3.561e-06  5.6899e-07
2.301e-06  7.719e-06   5.71e-07 
2.729e-06  1.4475e-05 6.1099e-07
3.636e-06  2.9217e-05 6.1901e-07
 6.1e-06   6.2638e-05 1.351e-06 
 9.45e-06  0.00012289 1.674e-06 
2.1552e-05 0.00027303 2.298e-06 
4.3468e-05 0.00049021 5.251e-06 
0.00012995 0.0008693  2.1542e-05
0.00031187  0.001089  4.4333e-05
0.0010387  0.0021593  0.00016535
 0.003756  0.0046594  0.00054881
 0.015004   0.010793  0.0021031 
 0.063103   0.02465   0.0068534 
 0.26945    0.062558   0.023848 
  1.1691    0.17201    0.083109 
  5.0574    0.47685    0.27839  
  26.491     1.6898     1.679   
  150.52     6.6698     4.4145  
  1062.5     19.547     14.176